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Introduction 

The solution of problems in most parts of applied mathematics and some areas of 
pure can more often than not be reduced to the problem solving some differential 
equations. Furthermore, many parts of pure maths were originally motivated by 
issues arising from differential equations, including large parts of algebra and much 
of analysis. From Mods, and even from school, you now know how to solve some 
differential equations. In this course, we shall learn to solve a much larger variety 
both explicitly and qualitatively, and we also develop some general theory. 

Throughout, we will use the following convenient abbreviations: we shall write 
DEs: for differential equations. 

ODEs: for ordinary DEs, i.e. differential equations with only ordinary derivatives. 
PDEs: for partial DEs, i.e. differential equations with partial derivatives. 
The course contains six topics, with a chapter devoted to each: 
The six chapters 

1. ODEs and Picard's Theorem (for existence of solutions). 

2. Boundary Value Problems (BVPs) for second-order linear ODEs. 

3. Plane autonomous systems. 

4. First-order quasi-linear PDEs: the method of characteristics. 

5. The classification of second-order linear PDEs. 
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6. Laplace and Fourier transforms (for ODEs and PDEs). 



The first three are on ODEs and the last three largely on PDEs. 
Books 

The main text is P J Collins Differential and Integral Equations, O.U.P. (2006), 
which can be used for the whole course (Chapters 1-7, 14, 15). 

Other good books which cover parts of the course include 

W E Boyce and R C DiPrima, Elementary Differential Equations and Boundary 
Value Problems, 7th edition, Wiley (2000). 

E Kreyszig, Advanced Engineering Mathematics, 8th Edition, Wiley (1999). 

G F Carrier and C E Pearson, Partial Differential Equations - Theory and Tech- 
nique, Academic (1988). 

J Ockendon, S Howison, A Lacey and A Movchan, Applied Partial Differential 
Equations, Oxford (1999) [a more advanced text]. 
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1 ODEs and Picard's Theorem 



1.1 What is an ODE? 

An ODE is an equation of the form 

G(x,y,y',y", y™) = 

regarded as an equation for y(x). Usually this can be solved for the highest deriva- 
tive and written in the form 

<1 " !l //" (.r) F(x,y,y', y^). 



dx r 



Then the order of the ODE is n, the order of the highest derivative which appears. 
Given an ODE, certain obvious questions arise. We could ask: 

• does it have solutions? Can we find them (explicitly or implicitly)? If not, 
can we at least say something about their qualitative behaviour? 

• given data e.g. the values y(a), y'(a), ... of y(x) and its first n — \ derivatives 
at some initial value a of x, does it have a solution? is it unique? Are other 
choices of data possible? 



We shall answer the first question in this section, and consider the second in the 
next section. 

For simplicity, we begin with a first-order ODE with data: 

V = f(x, y) with y(a) = b. (1.1) 

This is an initial value problem or IVP, since we are given y at an initial, or 
starting, value of x . You know how to solve a variety of equations like this, and 
perhaps expect a unique solution for a given initial value as you may not have 
encountered the following difficulty 

1.2 The warning example: 

Consider this IVP: 

y' = Zy 2lz ; 2/(0) = 0. (1.2) 
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So separate the variables (a mods technique): 

dy 



3y2/3 



dx, 



integrate and use the data to find 



(i) A solution y = x 3 ; 

(ii) But evidently there is another solution: by direct checking y = will do; 

(iii) Now we can see that there are infinitely many. Pick a, b with a < < b and 
define 

y = (x — a) 3 x < a 
= a < x <b 

= (x - b) 3 b < x 

The solution is not unique (in fact far from it, since we've found infinitely many 
solutions). So for uniqueness in the problem (1.1) we must impose conditions on /. 
We discover what these are in the process of proving that the solution exists. This 
will be the first existence theorem which you've encountered (mathematicians are 
unique in proving that certain things exist, and that certain things don't, which 
is usually easier). The proof is quite technical, certainly the most technical thing 
in the course. In particular, y ou need to remember from Mods the Weierstrass 
M-test for convergence of a series of functions. Once we get past the proof of this 
theorem, most of the rest of the course is methods of solution of particular classes 
of equation. 

To be precise then, we shall seek a solution in a rectangle R about the initial point 
(x, y) = (a, b), so suppose R is \x — a\ < h, \y — b\ < k as in figure 1.1: 
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Figure 1.1: The rectangle R 
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Our first assumption is that f(x,y) is continuous in R. Now integrate (1.1) 
from a to variable x: 

W)\l = [ X f(t,y(t))dt = y(x)-y(a) 

so rearranging 

y(x) = b+ f f(t,y(t))dt. (1.3) 

J a 

We have transformed the differential equation to an integral equation - the un- 
known y is given in terms of an integral rather than a differential. There is a 
general theory of these, but we only need to deal with the particular case of (1.3). 
The standard approach is to seek a solution by iteration or successive approxima- 
tion. 



1.3 Picard's method of successive approximation 

We start with an initial guess and then improve it. The guesses, or successive 
approximations or iterates, are labelled y n (x) starting with yo(x). Take 

y (x) = b 
y n +i(x) = b + f* f(t, y n (t))dt 

That is, we start with the simplest guess, that y equals its initial value, and at 
each stage substitute the current guess into the right-hand-side of (1.3) to get the 
next guess. We need to know if this process converges, and if it does whether 
it converges to a solution of the problem (1.3). Consider the differences between 
successive approximations: 

e (x) = b 
e n+1 (x) = y n+1 (x) - y n (x) 

and note that 

n 

Vn(x) = ^2e k (x). (1.6) 



We have written y n (x) as the sum of a series, (1.6) and we need this series to 
converge, so that we need the differences e n (x) to get small. With hindsight, this 
needs some assumptions and conditions on /, but here is the key idea: notice that 

e n+ i(x) = y n +i(x) - y n {x) 
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= [ X [f(t,y n (t))-f(t,y n ^(t))}dt 

J a 

and recall that the integral of the modulus of a function is greater than the modulus 
of the integral (because the function can be negative). Therefore 

I e n+ i(a;)| < 

Recall we want the e n to get smaller. This motivates a definition: 

1.4 The Lipschitz condition 

A function f(x,y) continuous on a rectangle R satisfies a Lipschitz condition with 
constant A if 3 real positive A such that 

\f(x, u) — f(x, v)\ < A\u — v\ for (x, u) G R, (x,v) e R. (1.8) 

This is a new condition on a function, stronger than being continuous but weaker 
than being differentiable. It turns out to be the right condition to make Picard's 
Theorem, which is the existence theorem we want, work. 

1.5 Picard's Theorem 

The ODE (1.1) 

y' = f(x,y) with y(a) = b. 
has a solution in the rectangle R : \x — a\ < h, \y — b\ < k provided: 

P(i): / is continuous in R, bounded by M (so \ f(x,y)\ < M) and Mh < k. 
P(ii): / satisfies a Lipschitz condition in R. 

Furthermore, this solution is unique. 

Proof 

We want the series in (1.6) to converge as n — > oo. We break the proof into a 
series of steps: 

(i) each y n (x) is continuous and its graph lies inside R; 

continuity follows from (1.4) by induction, continuity of / and properties of 
the integral; for the graph property, note that 

\y n+ i(x)-b\ < | / \f(t,y n (t))\dt\ 



f 

J a 



\f(t,y n ) - f{t,y n -i)\dt 



(1.7) 
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<\M I 

J a 



dt\ = Mix -a\< Mh < k 



where we've used P(i) and induction. Thus y n (x) doesn't get further from 
b than k, and so the graph of y n goes right across R (it 'comes out the side 
not the top or bottom'). 




Figure 1.2: successive iterates graphed in R. 

(ii) The Lipschitz condition P(ii) means that, for some positive A, 

\f(t,y n ) - /(t,y n -i)| < A \Vn - y n -i\ by (1.8). 

(iii) This enables us to prove that the e n (x) get smaller. From (1.7) and (1.8): 

|e n+ i(a;)|<| / \f(t, y n ) - f(t, y n -i)\dt\ 

J a 

<A\ \y n - y n -i\dt\ 



and so 



Now 



3n+l0*0l <^l / K(t)\dt\. 
J a 



(1.9) 



ei(x) = yx(x) - b = / f(t,b)dt. 

J a 

By P(i), / is bounded by M so that 

\ei(x)\ < | / \f(t,b)\dt\ < M\x-a\. 

J a 

Using this in (1.9) gives 

|e 2 (x)| < —^-\x-a\ , 



and then by induction, using (1.9) 

, A n ~ 1 M . .„ A n ~ x M ir> 

e n (x) < — x - a n < — h n . 

n\ n\ 

This is enough to make the series (1.6) converge: 

(iv) By the Weierstrass M-test, the series in (1.6) converges to a continuous 
function. 

This is immediate from (iii), since the series with n-th term A " ^f^" con- 
verges. Thus 

y n (x) ->■ y^x) 

with yoo(x) continuous. The iteration has a limit, which is a continuous 
function. Is this limit a solution of the original problem? 

(v) We can take the limit of (1.4) to find 

yM=b+ [ f(t, Voo (t))dt (1.10) 

J a 

(this needs justification: it uses a property of the integral and the conti- 
nuity of /) but the RHS in (1.10) is differentiable; therefore, so is y^x). 
Differentiate (1.10) to find 

y'oo( x ) = f(x,yoo(x)) 
and since also yoo(a) = b, it's a solution. Are there any more? 

(vi) Uniqueness: if y(x) and Y(x) are two solutions of the problem, consider their 
difference: 



i(x) = y(x) - Y(x) = [ X (f(t, y) - f(t, Y))dt 

J a 

so |e(x)| < | [ X \f(t,y)-f(t,Y)\dt\ 

J a 

<A\ f \y-Y\dt\ 



<A\ I \e(t)\dt\ 

J a 



where we have used the Lipschitz condition. Now e(x) is continuous on R; 
therefore, it is bounded say \e(x)\ < B so 



e(x)\ < A 



X 

Bdt 



AB\x - a\ 
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and inductively from this 

\ x _ a \ n A n h n 
I e(x)\ < BA n - — — < B — ; > as n ->■ oo and e(x) = 0. 

The difference is zero, so the solutions are the same i.e. unique (uniqueness 
proofs almost always go like this: assume there are two and make their 
difference vanish). 

This ends the proof of Picard's Theorem. Since the warning example doesn't have 
a unique solution, something goes wrong for it. As an exercise, show that the 
warning example fails the Lipschitz condition (in any neighbourhood of the initial 
point). 

1.6 Picard for first-order systems 

We've worked hard to prove existence and uniqueness for a single first-order ODE 
(f.f). Now we indicate briefly how the same technique can be used to solve a 
larger class of problems. Suppose then we have a system of first-order ODEs, (or 
a first-order system) in this case just a pair: 

y' = f(x,y,z) 

z' = g(x,y,z) 

with initial data 

y(a) = b; z(a) = c. 

We introduce a matrix notation (which is then clearly capable of dealing with 
systems of more than two ODEs): 




so that the system can be written as 

YL = E(x,Y); Y(a) = B 

We have written a system of any size in a way which is formally identical to (1.1). 
We can set about solving it by iteration, following (1.4) by setting 

Y n+l = B+ P F(t,Y n (t))dt. 
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Then this will converge as before if we have a Lipschitz condition 

\\F(t,u)-F(t,v)\\<A\\u-v\\ (1.11) 
with some definition of "size" for a vector, which could be 

2/i I + I2/2 j - 

Theorem: Picard for first-order systems 

The solution exists and is unique for F_ continuous and Lipschitz in a suitable 
region 

\y — b\ < ki, \z — c\ < k 2 , \x — a\ < h 
(and this extends to a system of n equations.) 

1.7 Picard for Higher Order ODEs 

With Picard extended to first-order systems, it is a small step to extend it to a 
single, higher order ODE. For simplicity, we consider just an I VP for a particular 
class of second-order ODEs, which are the concern of the next section where we 
learn to call them linear: 

y" + p(x)y' + q(x)y = 

with initial data 

y (a) = b y'(a) = c, 
and p, q continuous for \x — a\ < h. 

To reduce this to a first-order system, introduce z — y' and write 

y =z = f(x,y,z) 
z ' = - pz -qy = g(x,y,z) 

with data y(a) = b, z(a) = c. This is precisely in the form to which the previous 
section applies, and its easy to check that (1.11) is satisfied, so we get: 

Theorem: Picard for second-order linear ODEs 

With the assumptions we've made, the solution exists and is unique. 

Clearly this method can be extended to the IVP for an n-th order linear ODE. In 
particular, this justifies our belief that an n-th order ODE needs n pieces of data 
to fix a unique solution. 
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1.8 Things that can go wrong with ODEs 



We've proved existence and uniqueness of solution for a variety of situations. What 
kinds of thing can go wrong? 

(i) Non-uniqueness of solution, e.g. if Lipschitz fails. We have the warning 
example and it is easy to produce other examples. 

(ii) Blow up of solutions. Picard may give us a solution for a range of x but 
the solution may become singular (or more graphically blow up) if we try to 
extend x. As an example, the equation 

y' = y 2 

with data y(0) = b has solution y = -j^—; so y — > oo as x — > ^. 

b X 

For the second example, Picard will give existence in a rectangle but it may be 
small. E.g. with b = 1 and y' = f(x,y) = y 2 , choose h and k, fixing the rectangle. 
For the bound on /, M — (1 + k) 2 will do, and then to satisfy Mh < k we need 
h < (i+ k y ^ \- Therefore, for the proof of Picard with this initial point, we can't 
get past x = \. We could move to a new initial point and start again, but the 
explicit solution blows up at x — 1, so we can't get past that. 



1.9 Aside: The Interchange Trick 

This is a useful device, helpful with higher order ODEs. Suppose we are interested 
in the second-order (but not necessarily linear) ODE 

y" = f(x,y) 

with data 

y(0) = a, y'(0) = b. 

Integrate once, from to x: 



or 



y'(x)-b= f X f(t,y(t))dt 
Jo 



y\x)=b+ / f(t,y(t))dt 
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and then integrate again: 

px px pi 

y(x)-y(0)= bdt+ / f(s,y(s))dsdt 
Jo Jo Jo 



so 



px ft 

y(x) = a + bx + / / f(s,y(s))dsdt. 
Jo Jo 

Now in the double integral, interchange the order of integration: 



X ft 



f(s,y(s))dsdt= / f(s,y(s))dtds. 

JO Jo Js 




Figure 1.3: The region of integration for the interchange trick (shaded). 



Carry out the t-integration: 



r-x 

= / (x-s)f(s,y(s))ds 
Jo 



y(x) = a + bx + / (x — s)f(s,y(s))ds. 



so that 



This is the interchange trick, and it can be very useful. Here we've written the 
second-order ODE in terms of an integral equation (replacing two integrations by 
one) and we could imagine seeking a solution by iteration (though that is beyond 
this course). 

Having interchanged once we could imagine iterating that. Thus as an exercise 
show that, if 

r* ( x _ s )n-i 



y(x) 



o {n-l)\ 



-g{s)ds 
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then 

d n y , . 
= 9{x) 

and 

y(0) = i/(0)= ... = 2/ ("- 1 )(0)=0. 
This is the result of interchanging n — 1 times. 
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2 Boundary Value Problems for 2nd order linear 
ODEs 



In this chapter, we consider the problem of solving: 

y"(x) +p(x)y' + q(x)y = f(x), a<x<b 
y(a) = A, y(b) = B. 

This is a second-order ODE, so we need two bits of data, but here we are given y 
at the ends of an interval rather than y and y' at the same (initial) point: this is 
a boundary-value problem (BVP) not an IVP. We are seeking a solution through 
two specified points in the plane, namely (a, A) and (b,B). 

We want first to explain why this equation is linear. 

A second order linear operator is a map of functions 

y L[y] = P 2 (x)y" + P^y' + P (x)y. (2.1) 

We usually assume the coefficients Pi(x) are continuous, with P 2 (x) 7^ (i.e. not 
the zero function). This operator is linear in the sense that 

L[ a y 1 + f3y 2 }=aL[y l }+f3L[y 2 ] (2.2) 

for a, (3 G R (or C, but usually we stay with reals) and functions y±, y 2 - Thus its 
linear in the sense of being a linear transformation on a suitable vector space of 
differentiable functions. 

By a 2nd-order linear ODE we mean one of 

L[y\ := P 2 y" + P x y' + P y — : H, homogeneous 
L[y] := P 2 y" + P\y + P y = f : N, inhomogeneous 

for some given /. 
2.1 Solutions 

The following properties of solutions of H and iV are easily established: 
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(i) the solutions of H form a vector space (since if L[yi] = = L[y 2 ] then 
L[a yi + py 2 } =0). 

(ii) if y 1 and y 2 satisfy N then y 1 — y 2 satisfies H, so that the general solution 
of N may be written 

V = Vpi^ + Vcf, 

any solution of n general solution of h 

where ypj is called the particular integral and ycF the complementary function. 
Now we develop more of the vector space language. 

2.2 Linear independence of functions 

A pair of functions yi(x), y 2 (x) is linearly independent if the only linear combina- 
tion which vanishes (identically): 

ciyi(x) + c 2 y 2 (x) = 

has c\ = c 2 = 0. They are linearly dependent if such q, not both zero, can be 
found. If they are also differentiable then this would entail 

ciy[ + c 2 y 2 = 0, 

h y i)( Ci )=o 

\Vl 2/2/ V C 2/ 

so that the determinant of the matrix is zero. 

Define the Wronskian of a pair of function to be this determinant: 

W(x;y u y 2 ) = y 1 y 2 -y 2 y[. (2.3) 
From what we have just seen: 

2.3 Proposition 

If two functions are linearly dependent then their Wronskian vanishes. 
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The converse to this isn't obvious; consider the following pair of (once) differen- 
tiable functions: 






X 


< 


x 2 


X 


> o 


x 2 


X 


< 





X 


> 



then W = 0. However, if c\y\ + c 2 y 2 = then evaluation for positive and negative 
x shows that c\ = c 2 = 0, so that these functions are in fact linearly independent. 
To establish a partial converse note the following: 



2.4 Proposition 

If y 1 and y 2 satisfy 

H : P 2 y" + P iy ' + P y = 
then their Wronskian W satisfies. 

W ~ ~T 2 

The proof is an easy exercise. Solving for W, we get 

Pi(t) 



W = const x exp 



dt 



(2.4) 



In particular, provided P 2 is nowhere zero, if W = at one point, then W = 
at every point (since in this case the exponential can't vanish so this can only 
happen if the constant in front of (2.4) is zero). 



2.5 A basis of solutions of H 

We choose solutions yi and y 2 of H with 

yi(a) = 1, yi'(a) = 
y 2 (a) = 0, y 2 '(a) = 1. 

By the work in Section 1, these exist and are unique at least in a neighbourhood 
of x = a provided P 2 (a) 7^ 0. Also their Wronskian W(x) has W(a) = 1, so is 
nonzero in this neighbourhood of x — a, and they are linearly independent. Do 
they span the vector space of solutions? Suppose y(x) is any other solution and 

set 

Y(x) = y 1 {x)y(a) + y 2 (x)y'(a). 
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Then this is a solution with 



Y(a)=y(a); Y'(a) = y'(a) 



and so by uniqueness Y(x) = y(x). Therefore they do span the vector space of 
solutions, and therefore they are a basis. We can conclude: 

2.6 Proposition 

(i) The dimension of the space of solutions of if is 2. 

(ii) Any pair of solutions of H with W ^ is a basis. 



(Any linearly independent pair of solutions is also a basis, and so will be related 
to the basis of Section 2.5 by a nonsingular matrix. At a, their Wronskian will 
be the determinant of this matrix, and therefore nonzero. This gives our partial 
converse to Proposition 2.3: two solutions of a given fixed H are linearly dependent 
if and only if their Wronskian is zero. This fact can be proved directly as follows: 
suppose u and v are two solutions of H; if they are linearly dependent then we 
know already that their Wronskian is zero so now suppose for the converse that 
their Wronskian is zero; if u is the zero function then they are certainly linearly 
dependent so suppose that there is at least one value of x, say x — a, with u(a) ^ 0; 
pick jj, so that v(a) = /iu(a) then 



cancel u(a) to conclude that v'(a) = fj,u'(a); now define y(x) — v(x) — fj,u(x), then 
y(x) is a solution of H by linearity, while y(a) = = y'(a) by the choices made so 
far; thus by uniqueness of solution of H we conclude that y(x) = and therefore 
u and v are linearly dependent QED.) 

Exercise: generalise everything done so far to n-th order linear ODEs. 
2.7 Variation of parameters 

We now know a good deal about the solutions of H . 'Variation of parameters' is 
a method to use these to solve N . Recall the distinction: 



= W(a) = u{a)v\a) 



v(a)u'(a) = u(a)(v'(a) — fiu'(a)); 



L[y] := P 2 y" + P lV ' + P y 



: H 

f :N 
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and suppose that H is solved by y = ciyi(x) + c 2 y 2 [x) with linearly independent 
l/i, 1/2- We seek a solution of iV of the form 

y(x) = ci(x)yi(a;) + c 2 (x)y 2 (x), (2.5) 

i.e. we 'vary the parameters'. Using two functions to find one, we expect to be 
able to impose another condition on the q. 
Differentiate (2.5) to find 

y' = CiV'i + c 22/ 2 + c 'iVi + C 2V2 

and now impose 

c[yi + c' 2 y 2 = 0. (2.6) 

The justification for this is finally that it leads to an explicit formula for the 
solution of N. Differentiating again, 

y" = Cl y'( + c 2 y' 2 ' + d x y' x + c' 2 y' 2 

so that 

L[y] = P 2 {c iy 'i + c 2 y'i + c' iy [ + d 2 y' 2 ) 
+P 1 (c 1 y[ + c 2 y 2 ) 
+P (c 1 y 1 + c 2 y 2 ). 
But, since the yi satisfy H, this gives N as 

L[y] = P 2 (c' iy [ + J 2 y' 2 ) = f. (2.7) 

Solve (2.6) and (2.7) for c[ to find 

c i = —p^yi W = 2/13/2 - 2/22/i 



and then 



Integrate these to obtain 



/ fy 

c 2 " 



P 2 W 



c (x ) = _ r mv2{t) dt 
c n x ) J p 2 (t)w(t) aL 

^2\X) - J p 2{ t) W{t )UL 



(2.8) 



The freedom in the choice of lower limit in these integrals gives additive constants 
in Q and so, by (2.5), adds a solution of H to this solution of N. We shall see how 
to exploit this freedom after doing an example. 
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2.8 An example 

Consider the equation (an example of N) 

7T 

y" + y = tan a; for < x < — . 

The corresponding H is 

y" + y = 

for which we may choose two linearly-independent solutions as 

yi = sinx; y 2 = cosx. 
The Wronskian turns out to be 

W = y±y' 2 - y 2 y[ = -1 

so by (2.8) 

Ci — J tan t cost dt = ct\ — cos a; 

C2 = J — tan t sin t dt = a 2 — log(sec x + tan x) + sin x 
where the ctj are arbitrary real constants. Then the solution of iV is 

V = cm + c 2 y 2 
= — cos x log(sec x + tan x) + ct\ sin x + a 2 cos x . 

> . ' N v ' 

PI CF 

2.9 Fixing the boundary values: the Green's function 

Suppose we have the BVP 

Piy" + Piy' + P y = f (2.9) 

with data 

y(a)=0 = y(b). (2.10) 

We may solve this by Variation of Parameters, but choosing yi so that y\(a) = 0, 
and y 2 so that y 2 {b) = 0. (Assume this can be done; we shall return to this point.) 
So 

y(x) = ci(x)yi(x) + c 2 (x)y 2 (x) 
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with the Q as in (2.8). Then 

y(a) = ci(a)y 1 (a) + c 2 (a)y 2 (a) 

y(b) = ci(6)yi(6) + c 2 (b)y 2 (b). 

We want both these to vanish, and, with the choices made for y^ they will if we 
impose c 2 (a) = = Ci(fe). In this case, instead of (2.8) we have 

fb f(t)y 2 (t) 



Cl( v) = r ntmw dt 

.. i ..) _ r* MviVLjf < [ - n) 

^2\X) — J a p 2 (t)W{t) aL - 



Thus we fix the q completely. Now 

y(x) = ci(x)yi(a;) + c 2 (:r)w 2 (:r) 
[* f(t)Vi(t)V2(x) f b f(t)y 2 (t) yi (x) 

Ja P2(t)W(t) + ] x P 2 (t)W(t) 

which we can write concisely as 



dt 



y(x)= [ G(x,t)f(t)dt (2.12) 

J a 



where 

r yi{t)y2{x) a<t<x <b 

g(».«)= SO : ■ (2.i3) 

I p 2 (t)w(t) — — — 
We call G(a;,t) the Green's function for the BVP (2.9)-(2.10). 

Note that G(x, t) is defined and continuous on the square a < x < b, a < t < b in 
the (x,t)-plane, and, as a help to remembering the formula (2.13), observe that it 
vanishes on the sides of this square. (Although G(x,t) is continuous, dG(x,t)/dx 
has a jump across x — t of magnitude 1/P 2 (t).) 

2.10 An example 

Extending the example of section 2.8, consider the BVP 

7T 

y" + y = /(#) for < rr < — 

with data 

y(o) = o = y(|), 

and run through the method: 
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• Identify H, as before, as y" + y = 0. 

• Choose solution y 1 with yi(0) — so y\ — sinx will do. 

• Choose solution y 2 with y 2 (f ) = so y 2 = cosx will do. 

• Calculate 



W 



Vi V2 

y[ y'2 



-l. 



• Note P 2 = 1 
and so by (2.13) 



G(x,t) = 

The solution of the BVP is then, by (2.12) 



sin t cos x < t < x < | 
cos tsinx < x < t < | 



7T 

y(a;)= [ * G(x,t)f(t)dt, 
Jo 



with this G. 



(One can think of equation (2.12) metaphorically as follows: the BVP is L[y] = f 
so the solution must be y = L~ l [f] and this formula defines L^ 1 .) 



2.11 A warning 

This procedure will fail if W — 0, i.e. if y\ and y 2 and linearly dependent, i.e. 
proportional, which means that H has a (nontrivial) solution vanishing at both 
ends. Now think about uniqueness for the BVP: if the solution of iV is not unique 
then the difference between any two gives a solution of H satisfying both boundary 
conditions, which means vanishing at both ends. We conclude: 

• G is well defined iff 

• H+ boundary conditions has no nontrivial solution iff 

• N+ boundary conditions has a unique solution. 

If any one of these fails, they all do. If that happens, we must revert to Variation 
of Parameters as in section 2.7. 
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As an example, consider the problem 

y" + y = f(x), y(0) = = y(7r), (2.14) 

so H is y" + y = 0. 

We try to follow the method: choose y\ so that yi(0) = 0, e.g. yi = sin a;; but then 
yi(n) = too. H has a solution vanishing at both ends so that G is not defined. We 
revert to Variation of Parameters, choosing a second, linearly independent, solution 
y 2 of H, say y 2 = cosx. Then W — — 1, P 2 — 1 and y = ci(a;)yi(:r) + c 2 (x)y 2 (x) 
with Q as in (2.8). Now we try to impose the boundary conditions: 

y(0) = c 2 (0) therefore = and c 2 = — fit) sm tdt 

Jo 

by (2.8); but now 

y (n) = -c 2 (vr) 

and, for the other boundary condition we need this to vanish, which is the condition 

f{t)smtdt = 0, (2.15) 

a necessary condition on /, without which the original problem has no solution. 

We conclude that (2.14) has a solution iff / satisfies (2.15). If / does satisfy (2.15), 
then (2.14) has infinitely many solutions, since Ci(x) is undefined up to an additive 
constant. 

(Compare this with the situation of the matrix equation Ax = v where A is a 2 x 2 
matrix and c and v are two-component column vectors. If A is nondegenerate 
there exists a unique solution. If A is degenerate, it has a homogeneous solution 
Acq = and there is a column vector w such that w t A = 0. We can therefore 
only solve the equation when w l v = w 1 Ac = and then, if this is satisfied, given 
a solution c, it is not unique as c + Xcq is also a solution. The problem (2.14) 
can be reduced to this linear algebra problem by taking a particular integral Y 
to the differential equation and two independent solutions to the homogeneous 
equations y\ and y 2 and setting y = Y + C\y\ + c 2 y 2 for constants C\ and c 2 . Then 
the boundary conditions reduce to equation of the form Ac — v for the constants 




2.12 Example: Infinite limits 

Up to now we have considered BVPs on finite intervals. The extension to infinite 
intervals is straightforward but may need care with the limits at the ends of the 
interval. Here is an example: 
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Consider the following BVP: 



u" — u — f, for — oo < x < oo 

where / is bounded and continuous on the whole real line. Show that there is at 
most one solution which is bounded, and show that it can be written in the form 

/oo 
K(\x -t\)f(t)dt (2.16) 
-oo 

for a K which you should find. 



Solution: If u\ and u 2 are two bounded solutions of the BVP, consider u = wi— w 2 ; 
this is bounded and solves u" — u = 0; the general solution of this homogeneous 
equation is u = Ae x + Be~ x but this solution can only be bounded for all x if 
A = B = 0. That proves uniqueness. 

For the explicit solution, construct the Green's function choosing y 1 — > as x — > 
—00 and i/2 — > as x — > +00; so yi = e x and y 2 = e~ x will do. Then 

W = -2, P 2 = 1 
1 



so 



Cl = "2 



C2 = "2 



and 



Define 



u = 



fitydt 

CO 



e*- x f(t)dt 



e x ~ l f{t)dt. 



K = — -e-l*-*" 
2 

then this takes the form of (2.16). This is a solution, but is it bounded? Suppose 
\f \ < M then 



\u(x)\ < 



\f(t)\dt 



< M. r e -\*-t\ dt 



00 
00 



< M 



e~ l dt 



M 



so u is bounded, and is therefore the unique bounded solution, by the first part. 
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3 Plane autonomous systems of ODEs 



The definition: a plane autonomous system of ODEs is a pair of ODEs of the form; 

ft T 

Here "autonomous" means there is no t-dependence in X or 7, and "plane" means 
there are just two equations, so we can draw pictures in the (x, y) - plane, which 
will then be called the phase plane. Given initial values x(0) = a, y(0) = b, 
we expect to find a unique solution (this needs X and Y Lipschitz, which we 
therefore assume) which will define a trajectory or phase path in the phase plane. 
It is convenient, though not necessary, to think of t as time, and the trajectory as 
tracing out the path in time of a moving particle. Then we can put an arrow on 
the trajectory giving the direction of increasing t. 



3.1 Critical points and closed trajectories 

A critical point is a point in the phase plane where X = Y = 0. 
A critical point is a particular (very special) trajectory (for which the particle 
doesn't move). Note that trajectories can only cross at a critical point (since 
otherwise the particle has more than one velocity there). 

There may be trajectories in the phase plane which are closed i.e. which return to 
the same point. Provided they dont contain a critical point, these correspond to 
periodic solutions of (3.1) as may be seen as follows: 

suppose that for some finite value p of t, (x(p),y(p)) = (x(0), y(0)), while (x(t),y(t)) ^ 
(x(0), y(0)) for < t < p. Define x(t) = x(t + p), y(t) = y(t + p). Then 

dx 

— = X(x(t + p), y(t + p))=X(x(t), y(t)) 

and 

f f = Y(x(t), y(t))- 

So (x(t),y(t)) is another solution with x(0) = x(p) = x(0); y(0) = y(p) = y(0). 
Now by uniqueness of solution (given Lipschitz again). 

x(t + p) = x(t) = x{t) 
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y(t + p)=y(t)=y(t), 



but this is now true for all t, so a closed trajectory corresponds to a periodic 
solution of (3.1) with period p. (It would be possible to have infinite p but then 
the trajectory has a critical point in its infinite past and the same critical point in 
its infinite future.) 

3.2 An example 

Consider the harmonic oscillator equation 

x = —co 2 x. 

Turn this into a plane autonomous system by introducing y as follows: 

x = y =X(x,y) 1 . 32 v 
so y = —u 2 x =Y(x,y). J 

(Clearly this trick often works for second-order ODEs arising from Newton's equa- 
tions.) The only critical point is (0,0), but note that 

^-(u 2 x 2 + y 2 ) = 2u 2 xx + 2yy = 

so oo 2 x 2 + y 2 = constant, (which, from Mods mechanics, we know to be propor- 
tional to the total energy). For a given value of the constant this is the equation of 
an ellipse, so we can draw all the trajectories in the phase plane as a set of nested 
(concentric) ellipses: 

yk 




Figure 3.1: The phase portrait for the harmonic oscillator; to put the arrows on the trajec- 
tories, notice that x > if y > 0. 
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The picture in the phase plane is called the phase portrait and from that we see 
that all trajectories are closed, so all solutions are periodic (as we already know, 
from Mods). 

We want to learn how to sketch the trajectories in the phase plane in general. We 
start by classifying critical points, so suppose P = (a,b) is a c.p. for (3.1), so 

X(a,b) = = Y(a,b). (3.3) 

Now x = a, y = b is a solution of (3.1). We analyse its stability by setting 

x = a + C(t); y = b + rj(t) 

where ( and 77 are thought of as small. From (3.1), 

x = C = X(a + C, b + rj) = X(a, b) + (X x \ p + r]Xy\ p + h.o. 

y = 77 = Y(a, b) + (Y x \ p + r]Yy\ p + h.o. 

where 'h.o.' means quadratic and higher order terms in ( and 77. Now use (3.3) 
and neglect higher order terms to find 



(c d) (J 

^ 1 1 ^Cr I p -^y I p 



(3.4) 



Call this (constant) matrix M and set Z_(t) = then (3.4) becomes 

Z = KZ. (3.5) 

We can solve (3.5) with eigen-vectors and eigen-values as follows: Z_ e xt is a solu- 
tion, with constant vector Z_ and constant scalar A if 

\z = mz , 

i.e. Z_ is an eigen-vector of M with eigen-value A. We are considering just 2x2- 
matrices, with eigen-values say Ai and A2 so the general solution if Ai 7^ A2 is 

Z(t) = ciZ ie Al * + c 2 Z 2 e X2 \ (3.6) 

for constant q. Recall Ai, A 2 may be real, in which case the c, and the are real, 
or a complex conjugate pair, in which case the q and the Z_ t are too. 

If Ai = A 2 = A G K say, we need to take more care. The Cayley-Hamilton Theorem 
implies that (M — A/) 2 = since the characteristic polynomial is c M (x) = (x — \) 2 , 
so either M — XI = or M — XI ^ 0. We have a dichotomy: 
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(M - XI)Z — (M — \I) 2 Z_ X = 0. 
One now checks that the solution of (3.5) is 

(ci£i + (c + C!t)Z )e A '. 



(3.7) 



Now we can use (3.6) and (3.7) to classify critical points. 
3.3 Classification of critical points 

We shall assume that neither eigenvalue of the matrix M is zero, which is the 
requirement that the critical point be non- degenerate. A proper discussion of this 
point would take us outside the course but roughly speaking if a critical point is 
degenerate then we need to keep more terms in the Taylor expansion leading to 
(3.4), and the problem is much harder. 

Case 1. A 2 > Ai > (both real of course) 

From (3.6), as t — > — oo, Z_(t) — > 0, while as t — >■ +oo, Z_(t) ~ a large multiple of 
Z 2 , unless C2 = when Z_{t) ~ a large multiple of 




Figure 3.2: An unstable node. 
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These trajectories converge on the critical point into the past, but go off to infinity 
in the future. A critical point with these properties is called an unstable node 
or an unstable improper node. 



Case 2: Ai < A 2 < (both real) 

This is as above but with t — ?► — t and the roles of Z 1 , Z_ 2 switched. The trajectories 
converge on the critical point into the future and come in from infinity in the past. 




Figure 3.3: A stable node. 
This is a stable node or stable improper node. 

Case 3: Ai < < A 2 (both real) 

If ci = then Z_(t) — > oo along Z_ 2 as t — > oo 

— > as t — > —oo. 

If C2 = then Z_(t) — > along Z_ x as t — > oo 

— ?■ oo as t — > — oo. 

If ci, c 2 7^ then Z_(t) — > oo along Z_ 2 as t — > oo along Z_ x as t — )■ — oo. 

Most trajectories come in approximately parallel to ±Z 1 and go out asymptoting 
to ±Z_ 2 . 
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Figure 3.4: A saddle. 



This is a saddle (to motivate the name, think of the trajectories as contour lines 
on a map; then two opposite directions from the critical point are uphill and the 
two orthogonal directions are downhill). 

Case (i) of equal roots is a proper node, stable if A < and unstable if A > 0. 
Case (ii) of equal roots is again an improper node but the phase portrait is different. 




Figure 3.5: Unstable proper node case (i) and unstable improper node case (ii) 

If the eigen- values are a complex conjugate pair we may write 

Ai = /i — iv, A 2 = n + iv /i,i/Gl, 

and the classification continues in terms of /i and v. 
Case 4: /x = 
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so Ai = —iv and A^ = — v 2 < 0; in terms of the matrix M of (3.4), A + D = but 
AD - BC > so BC < 0. Equation (3.4) becomes _ 



A 
C 



B 
-A 



(3.8) 



As an exercise, show that now —C( 2 + 2A(i] + Brf is constant in time. We know 
that B,C have opposite signs with (—BC) > A 2 so this is the equation of an 
ellipse. 




Figure 3.6: A clockwise centre (B > 0) 

This case is called a centre. The sense of the trajectories, clockwise or anticlock- 
wise, depends on the sign of B; as can be seen from (3.8), B > is clockwise (take 
C = and 7] > 0, then ( = Br] > 0). 

Case 5: ix ^ 

So, in (3.6), we must have Z_ x = Z_ 2 and c\ = C2 and 

Z(t) = [c x Z x e- ivt + Cil^'] , 

which is just like case 4, but with the extra factor e M *, which is monotonic in time. 
We have another dichotomy: 

(i) ix > then \Z_(t)\ — > oo as t — > oo so the trajectory spirals out, into the 
future. This is called an unstable spiral. 
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Figure 3.7: An unstable spiral; reverse the arrows for a stable spiral 



(ii) fj, < this is the previous with time reversed so it spirals in, and is called a 
stable spiral. 

In case 5, as in case 4, the sense of the spiral is dictated by the sign of B. Impor- 
tant observation: if A + D > then we have one of the cases 1, 3 or 5(i), all of 
which are unstable (but if A + D < the critical point can be stable or unstable). 



3.4 An example 

Find and classify the critical points for the system 



solution: for the critical points, from X = deduce x — y, therefore from Y = 
deduce x 2 = 1, and we have either (1, 1) or (—1, —1). 

For the classification, calculate 



x = x — y = X(x, y) 



y = 1 - X y = Y(x,y) 




and evaluate at the cps: 



at (1,1) : M 



1 



-1 
-1 



A 2 -2 = 0: A = ±V2 



this is a saddle. The corresponding 



eigenvectors are: 
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A2 = V2 Z_ 2 = ^ ^ y/^J direction ou ^ 



at (-1, 1) : M = Q : A 2 - 2A + 2 = : A = 1 ± %. 

this is an unstable spiral; B < 0, so its described anticlockwise. 




Figure 3.8: The phase portrait 



3.5 Further example: the damped pendulum 

Another example from mechanics: a simple plane pendulum with a damping force 
proportional to the angular velocity. We shall use the analysis of plane autonomous 
systems to understand the motion. 

Take 9 to be the angle with the downward vertical, then Newton's equation is 

ml9 = —mg sin 9 — mkl9, 
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where m is the mass of the bob, / is the length of the string, g is the acceleration 
due to gravity and A; is a (real, positive) constant determining the friction. We 
cast this as a plane autonomous system in the usual way: set x = 9 and y — x — 9 

so 

x = y 

y = —— situ — ky 

For simplicity below, we'll also assume that k 2 < y, so that the damping isn't too 
large. 

To sketch the phase portrait, we first find and classify the critical points. The 
critical points satisfy y = = sinrr, so are located at (x,y) = (Nn,0). Then 

M=( ,° \ 

\— |cosx — k 

The classification depends on whether N is even or odd: 

for x = 2mr M = 
which gives a stable spiral (clockwise); 

for x = (2n+ l)n M = 

which gives a saddle. 
We now have enough information to sketch the phase portrait (note that x is 
positive or negative according as y is). 




Figure 3.9: The phase portrait of the clamped pendulum 
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3.6 An important example: The Lotka— Volterra equations 

This is a simplified mathematical model of a predator-prey system. Think of vari- 
ables x standing for the population of prey, and y for the population of predators, 
both functions of t for time. As time passes, x increases as the prey breed, but 
decreases as the predators predate; likewise y increases by predation but decreases 
if too many predators compete. We assume that x and y are governed by the 
following plane autonomous system: 



x = ax — jxy (3.9) 

y = -fiy + Sxy, 

where a, (3, 7, 5 are positive real constants. Because of the interpretation as 
populations, we only care about x > 0, y > but we shall consider the whole 
plane for simplicity. Again, the aim is to use the analysis of plane autonomous 
systems to lead us to the phase portrait and an understanding of the dynamics. 

For the critical points first, set 

X := x(a — 7?/) = 

Y :=y(-/3 + 5x) = 0. 
There are two solutions, (0,0) and (f, ^). For the matrix: 

m = ( Xx x y\ - ( a ~ 1V ~ 1X \ 
\Y X Y y J \ 5y -(3 + Sx) 

so first 

at (0,0): M= (jj 

which gives a saddle, where, it is easy to see, the out-directon is the x-axis and 
the in-direction is the y-axis. Next 



<T 7/ ' Vf 



.Pi 



ai (4, -V M= (i J ) : A 2 + n 



which gives a centre, described anticlockwise since B < 0. 

We have found and classified the critical points. Before sketching the phase por- 
trait, it is worth noting, from (3.9), that if x = then x — and if y = then 
y — 0. Thus the axes are particular trajectories, and trajectories can only cross at 
critical points (as noted before). 
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Figure 3.10: The phase portrait for the Lotka-Volterra system 

Therefore any trajectory which is ever in the first quadrant is confined to the first 
quadrant, and no trajectory can enter the first quadrant from outside. Since there 
is a centre in the first quadrant, it looks as though all trajectories in the first 
quadrant are periodic. This is true, and can be seen by the following argument: 
form the ratio 

y = y(-fi + ax) = dy_ 
x x(a — 72/) dx 

and separate 

{a - iy) dy - tl±M dx = o ; 

y x 

now integrate 

(3 log x — 5x + a log y — jy = const 

or 

xPy a e- Sx e-™ = C, (3.10) 

for a constant C, which is necessarily positive for a trajectory which is ever in the 
first quadrant. For different values of C, (3.10) is the equation of the trajectory 
or equivalently the trajectories are the level sets or contours of the function on 
the left in (3.10). This function has a single maximum, vanishes on the axes and 
goes to zero at infinity. Therefore its contours are all closed curves and so all the 
trajectories are closed and all the solutions of (3.9) are periodic. 

This useful technique can be applied to other examples. 
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3.7 Another important example: limit cycles 



Consider the plane autonomous system: 



x = (1 - {x 2 + y 2 ) 2 )x - y 
y = (l- (x 2 + y 2 )^)y + x. 



Without the first term in each, this is the harmonic oscillator. 
Put x 2 + y 2 = r 2 then 

X = x(l — r) — y 
Y = y(l — r) + x 

and one sees that only critical point is (0, 0). One can go through the classification 
for this, to find that it is an unstable spiral (exercise!) or discover this fact from 
the following. 

We shall transform to polar coordinates. The simplest way to do this is as follows: 
first 



rr = xx + yy = x[x(l — r) — y] + y[y(l — r) + x] 
= r 2 (l -r) 



or 



f = r(l — r). 



Then, with 



y = rsin#, 



we find 



y = r sin 6 + r cos 66 = y(l — r) + x, 



which gives 6, so the system becomes 



6 = 1 



f = r(l — r). 



Unlike the system in its previous form, we can solve this. First 



6 = t + const, 



and then 




so 



log 



r 



1-r 



— t+ const 
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Solve for r and change the constant: 

1 1 

'' " 1 + Be-* ~ 1 + (i - l)e-* 

where r(0) = r . 

Note that as t — > oo, r — > 1, while as t — > — oo either r — )> if r < 1 or r — )> oo 
at some finite t if r > 1. 

Now it is clear that the origin is an unstable spiral, and that the trajectories spiral 
out of it anticlockwise. We can also see that r = 1 is a closed trajectory and 
that all other trajectories (except the fixed point at the origin) tend to it; we call 
such a closed trajectory a limit cycle. It is stable because the other trajectories 
converge on it. (For an example of an unstable limit cycle we could consider the 
same system but with t changed to —t.) 




Figure 3.11: Phase portrait with a limit cycle 



Another system with a limit cycle arises from the Van der Pol equation: 

x + e(x 2 — l)x + x = 
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where e is a positive real constant. If e = this is the harmonic oscillator again. 
If e 7^ then the usual trick produces a plane autonomous system: 

x = y 

y = —e(x 2 — l)y — x. 
The only critical point is (0, 0) and its an unstable spiral for e > (exercise!). 

Claim: Its beyond us to show this, but this system has a unique limit cycle, which 
is stable. There are some good illustrations for this in e.g. Boyce and di Prima 
(pp 496-500 of the 5th edition). 



3.8 Conservative and Hamiltonian systems 

Newtons law in one dimension with a conservative force 

dV{x) 



x = 



dx 



has a conserved energy H = x 2 /2 + V(x). This can be expressed as the plane 
autonomous system 

dVix) 

x = y, y = — ~ — , 
ox 

with conserved quantity H = y 2 /2 + V(x). More generally this is an example of a 
Hamiltonian system 

dH dH 
= ~k i y = 



dy ' ' dx ' 
where H := H(x,y) is an arbitrary function. Given an arbitrary system 

x = X(x,y), y = Y(x,y) 

it is both necessary and sufficient that 

d dX dY _ 
dx dy 

for the system to be Hamiltonian (necessity follows from the fact that partial 
derivatives commute, and sufficiency follows from the fact from mods that this 
equation is the integrability condition that allows us to integrate X — H y and 
Y = -H x for H). 
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A critical point for such a plane autonomous system is necessarily a critical point 
for the function H (i.e., its first derivatives vanish). At such a critical point, the 
matrix of the linearization is 

M = (-H -H ) 

and in particular it is trace free. The classification of critical for the equation 
reduces to that of critical points for H . In particular, nondegenerate critical points 
are either maxima, minima or saddle points for the function H as learnt in mods. 
The saddle is case 3 above, but the trace-free condition implies that Ai = — A 2 . 
Maxima and mininima correspond to centres. 

A key property is that H is constant along the flows (exercise) and so the flows can 
be thought of as contours of the graph of H over the plane. This is very useful for 
drawing the flows in the phase plane, and we can, for example, immediatly infer 
that flows in the neighbourhood of a centre are closed as they are contours near a 
maximum or minimum of H. 

Clearly these are very special and well-behaved systems. 



3.9 The Bendixson— Dulac Theorem 

It's important to be able to detect periodic solutions, but it can be tricky. We end 
this section with a discussion of a test that can rule them out. 

Consider the system x = X(x,y), y = Y(x,y). If there exists a function ip(x,y) 
with 

in a simply connected region R then there can be no nontrivial closed trajectories 
lying entirely in R. 



Proof. By nontrivial, I mean I want such a trajectory to have an inside i.e. it 
isn't just a fixed point. So suppose C is a closed trajectory lying entirely in R and 
let D be a disc entirely in R whose boundary is C. Consider the integral 




pdxdy 



D 




D 



dxdy 



= <j> —ipYdx + ipXdy 
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<j> —ip (—ydx + xdy) 



but on C, dx = xdt, dy = ydt so this is zero, which contradicts positivity of p, so 
there can be no such C. 



3.9.1 Corollary. 

If 

dX dY 
dx dy 

has fixed sign in a simply connected region R, then there are no nontrivial closed 
trajectories lying entirely in R. 

This is just the previous but with ip const — in an example, always try this first! 



3.10 Examples 

(i) the damped pendulum (section 3.5) 

x = y 

9 . , 
y = —- sim — ky 

has no periodic solutions. 
Here 

OX dY , 

jr + jr = - k <0 ' 

ox ay 

now use the corollary. 

(ii) 

x + f{x)x + x = 

has no periodic solutions in a simply connected region where / has a fixed 
sign. 

By the usual trick we get the system 

x = y 

y = -yf(x) - x 
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then 

dX dY 

+ = - fix) 

and we use the corollary, 
(iii) The system 

x = y 



y = -x - y + x 2 + y 2 



has no periodic solutions. 



The corollary doesn't help so try the general case: 

p := (ipX) x + { V Y) y = (p(-l + 2y) + X(p x + Y<p v . 

Now guess: <p y — then p = <f(—l + 2y) + yip x so take <p = —e~ 2x and 
p = 2e~ 2x > and we are done. 

Clearly, if ip isn't constant we need luck or a hint! 



41 



4 First-order quasi-linear PDEs: the method of 
characteristics 



4.1 The problem 

In this chapter, we are interested in a PDE of the following form: 

dz dz 
P(x, y,z)— + Q(x, y,z)— = R(x, y, z) (4.1) 

This is clearly first-order, and is quasi-linear by being linear in the highest order 
partial. Note that P, Q are allowed to depend on z, so that neither the whole 
equation nor even the left-hand-side are actually linear in z. We shall assume that 
P, Q, R are continuous and Lipschitz in the region of interest, yet to be defined. 

There are two problems: find the (or a) general solution or find a unique solution 
given data and determine its domain of definition. This is the region in the (x, y)- 
plane in which the solution is uniquely determined by the data. It turns out to 
depend on both the equation and the data. 

The solution of (4.1) will be a function 

z = f(x,y) 

but can be thought of as the surface defined by this equation, or equivalently 
defined by the equation 

E(x,y,z):=z-f(x,y) = 0. (4.2) 

We shall refer to this as the solution surface and call it E. The method of solution 
of the equation will be to generate E. 





Figure 4.1: The solution surface 
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A normal to the solution surface is defined by 



(this is a fact from Mods for the single-subject mathematicians; Maths and Comp 
students should ask their tutors for enlightenment) so consider the vector t = 
(P,Q,R). Then 

t .n = -P°L-Q°L + R , 
ox ay 

which vanishes by (4.1), so t is tangent to the surface S. 



4.2 The big idea: characteristics 

We look for a curve T whose tangent is t. If F — (x(s),y(s), z(s)) in terms of a 
parameter s (not necessarily path- length) this means 

% = P(x,y,z), ) 

% = Q(x,y,z), \ (4.3) 
% = R(x,y,z). J 

These are the characteristics equations and the curve T is a characteristic curve 
or just a characteristic. Given a characteristic (x(s),y(s),z(s)), call the curve 
(x(s),y(s),0), which lies below it in the (x, y)-plane, the characteristic trace or 
characteristic projection. 

The next result shows that characteristics exist, and gives the crucial property of 
them: 

4.3 Proposition 

(a) Through each point of space there is a unique characteristic. 

(b) Given a point p G T, the characteristic through p lies entirely on S. 

Proof 

(a) This, with our assumptions on P, Q, R, follows from existence and uniqueness 
for systems. 
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(b) (This is fiddly.) Consider the function E as a function along T; if it stays 
zero then T lies on the surface E = 0. So 

Z(s)=z(s)-f(x(s),y(s)) 

whence 

— = R(x(s),y(s),z(s)) - —P(x(s),y(s),z(s)) 



—JjjQ(x{s),y(s),z(s)) 



using (4.3), while 



df df 
R(x, y, f) - ^P(x, y, f) - -±Q(x, y, f) = 0. 

using (4.1). Subtract these 

<9E df 
— = R(x, y, z) - R(x, y, f) + (P(x, y, z) - P(x, y, f)) — 

df 

-(Q(x, y, z) - Q(x, y, /))^-- 

Now put z — f + E, then the RHS here is a function F(E(s), s) for some F 
with F(0, s) = 0. 

We have an ODE for E(s), and the RHS is Lifschitz (check!), so it has a 
unique solution, but E(s) = is a solution, so it's the only one. Therefore, 
E = all along T, so T lies on the solution surface. Q.E.D. 

Thus the solution surface E is ruled or generated by a collection of characteristics. 

4.3.1 Examples of characteristics 

We need to gain proficiency in calculating characteristics. 

(a) Calculate the characteristics for the PDE: 

dz dz 
ox ay 

From (4.3) write down the characteristic equations and solve them: 

, s 

— = P = x; x = Ae s 
ds 
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^- = Q = y; y = Be s 
as 

<^L = R = Z . z = Ce s 
as 

with A, B, C constants (trivial to solve), 
(b) Calculate the characteristics for the PDE: 

, dz dz 

This is a little trickier: 

dx . . C _, 

— = x - z; x = Ae H e 

as 2 

dy „ 
/ = 1; y = B + s 
ds 

^ = —z\ z = Ce~ s . 
ds 

with A, B, C constants. To solve this system, pass over the first, solve the 
second and third, then come back to the first. (I am adopting a conven- 
tion to introduce the constants A, B, C in the first, second and third of the 
characteristic equations respectively.) 

Solving the characteristic equations needs experience and luck; there isn't a general 
algorithm. 



4.3.2 The Cauchy problem 

The Cauchy problem will be to identify initial data that determines a unique 
solution on some domain and to identify the domain on which the solution is 
uniquely determined. 

Suppose we are given the solution z of (4.1) along a curve 70 (the data curve) in 
the (x, y)-plane. This produces a curve 7 in space: 
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y(t) = (x(t),y(t),z(t)) 




7o (i) = (x(t),y(i),0)) 
Figure 4.2: Geometry of the Cauchy problem. 



We introduce a parameter t along 7 so it is (x(t),y(t), z(t)) while 70 is the pro- 
jection (#(£), y(t), 0). Then, to solve (4.1), we construct the solution surface S by 
taking the characteristics through the points of 7. Thus the method of solution, 
the method of characteristics, is 

(i) Parametrise 7 as (x(t),y(t), z(t)). 

(ii) Solve 

ds 

7T = « 
as 

as 

for (x(s,t),y(s,t), z(s,t)) with data rr(0,£) = 2/(0, t) = y(t); z(0,t) = 

*(*)■ 

then, knowing (x(s, i), y(s, i), z(s, t)), we have found E parametrically i.e. as a 
function of s and t. 

We would like the solution explicitly, that is z in terms of x and y, a question to 
explore below, and there is a restriction on the data for the method to work, also 
to be found later. 
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4.4 Examples 

(a) Solve 



dz dz 



with z — 1 on x — y for < x < \. 



A 


7 












- — 7o 




1/2 ^ 



Figure 4.3: The data curve for this problem 



We introduce a parameter t for the data, say j(t) = (t, t, 1), for < t < 1/2, 
and then solve the characteristic equations (done in section 4.3.1) with this 
as data at s = 

C C 
x = Ae s + —e~ s ; x(0, t) = A + — = t 

y = B + s; y (0,t) = B = t 
z = Ce- s ; z(0,t) = C=l 
So, C — \, B — t, A — t — \ and the parametric form of the solution is 



x = (t- \)e s + \e- s 
y = t + s 



(4.4) 



z = e 



for < * < \. 



(b) Solve 



dz dz 
dx dy 
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with z = sinx on y = 0. 

So we can take ^(t) = (i, 0, sint), and solve the characteristic equations: 
= P = Z] x = A-2Ce—* s ; A - 2C = t 

as 



dy_ 

dx 



Q = l; y = B + s; B = 



- = R=--z; z = Ce 



C — sin t 



where the third column is the value at s = 0. So the parametric form of the 
solution is 



x = (t + 2sint) — 2sinte 2 s 
y = s 

z — suit e 2* 



(4.5) 



4.5 Domain of definition 



Where is the solution determined uniquely by the data? This is the domain of 
definition and will be the region in the (x, y)-plane where the solution surface is a 
graph of the function z = f(x,y). The solution surface may have edges and may 
have singularities. We want z as a function of x and y, so we need to be able to 
eliminate s and t in favour of x and y, at least in principle. For this, recall from 
Mods the definition of the Jacobian: 



J 



d(x,y) 



d(s,t) 



det 



%s Us 

xt Vt 



(4.6) 



The significance of J, again from mods, is that it relates the area elements: 

drr dy = Jds dt 

(this may not be known to Maths and Comp students; ask your tutor for enlight- 
enment). Now a necessary condition to be able to eliminate s and t in favour of x 
and y is that J be finite and non-zero, because we need non-zero infinitesimal areas 
on the plane and on the surface to correspond (think about this geometrically); in 
simply-connected regions of the (x, |/)-plane this condition on J is also sufficient - 
that's harder to prove so we'll just believe it for now. 

Let us determine the domain of definition for the example of section 4.4(a): 

The solution surface has edges given by the characteristics through the ends of 7, 
which are at t — and t — 1/2. 
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At t — 0, the characteristic trace is 

x = — sinhs; y = s so x = — sinh?/; 



at t — \ it's 



1 -s 1 1 ] 




x = — sinh y 

Figure 4.4: The domain of definition for this problem 



so the surface lies above the region 



— sinh y < x < —e 2 y 
y - - 2 



and these curves don't meet (check!). Next 

J _ I. ,4- / X S Us 



det 



det 

" x t y t 
(t - |)e s j - 



1 
1 



2> 2 

s 



(4.7) 



a - \y - \e~> 

but t — | < 0, so this is never zero (or oo), and the surface has no singularities. 
Therefore (4.7) gives the domain of definition. 
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Now look at the example of section 4.4(b); there are no edges but 



J 



det 



* 1 
x t 



Ft 



= |l + 2cost(l -e-2 s )\, 
where * is unimportant for the determinant. 

This is non-zero for s = i.e. initially, but it can vanish. Consider the range 



log - < s < log4 



so 



, 3 1 , n 

lo g2 > > ~ log2 

3 _i 8 1 
- > e 2 > - 
2 2 

-1 < 2(1 -e"5 s ) < 1 

then in this range 1 + 2 cost(l — e~^ s ) > 0. Note y = s here so for 

3 

-2 log - < y < 21og2 

the surface E is well-defined, but outside this range £ "curls up" (J vanishes and 
the surface becomes vertical before overlapping itself, ceasing to be a graph). 



4.6 Cauchy data 

We could calculate J at s = i.e. at 7. There it is 



det 



X s Vs 

x t y t 



det 



P Q 

x y 



using a dot for d/dt at 70 = (x(t),y(t),0). So on 70 we require 

P(x, y, z)y - Q(x, y, z)x ^ 0. 



(4.8) 



This is a necessary condition on the data for the solution to exist; if it holds, 
call the data Cauchy data. (Geometrically, the condition is that the characteristic 
traces at 70 are not tangent to it.) 
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4.7 Multivalued solutions and blow-up 



To get a better grip on what happens when J is or oo, we shall do another 
example, which is basic for this phenomenon: solve 



(4.9) 



z = f(y) on x = 0. 

(This is equivalent to the evolution equation z t + zz x = for a nonlinear wave with 
height z = z(x, t) whose wave velocity is determined by the height z, as Zt + cz x = 
is the equation for a wave with constant velocity c, solution z — z(x — ct).) 

We start with the characteristic equations: 

Ox 

— = p = l- x = s + A; X = A = 

OS 

^- = Q = z; y = B + Cs; y = B = t 

OS 

Oz 

R = ; Z = C; z = C = f(t) 

OS 

where, as usual, the third column gives the values at s = 0. Note that z is constant 
along any characteristic. 

The solution is 

x = s 

y = t + sf(t) 
* = f^) 

which we can rearrange as y — xz = t; z — f(t) or, by eliminating t, 



z = f(y- xz )- 

This is an important formula, giving the implicit solution of (4.9), and a lot of 
information can be obtained from it. 

Note that 

! + «/'(*) I • 

At s = 0, J = 1 7^ 0, so by (4.8) this is Cauchy data, but if /' < for any t, then 
J = at s = -^ij > 0. If we think of x as time, so that the data was given at time 
zero, then a problem has arisen at a finite time after the start. 



J 



det[ Xs Vs 

x t Vt 
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What is the nature of this problem? Calculate the partial derivative by the chain 
rule: 

dz dz dx dz dy ^z , 

so that 

dz f{t) 
dy 1 + sf'(t) 

and this goes to oo as s — > —jr^j, i- e - || ~^ 00 where J — > 0. Hence also 

dz dz 
dx dy 

This is known as a "finite-time blow up" of the partial derivatives (we don't know 
that z itself is singular). 

For more details, we look at a special case: consider the following / 

f(y) =0 y < 
-y 0<y<l 
-1 y>i 

We draw the characteristic traces, and recall that z is constant along the charac- 
teristic (figure 4.5). 
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Figure 4.5: The characteristic traces; note the convergence at x = 1 
From this we can sketch the graph of z(x,y) at x — 1/2, 1 and 3/2 (figure 4.6). 




Figure 4.6: The graph of z(x, y) at x = 1/2, 1 and 3/2 



From the figures, for x > 1, z is multi-valued, and the solution is not uniquely 
determined. 

In the context of this example, note that: 

1. we could do a smooth version with e.g. f(y) = —1 — tanhy and get a 
qualitatively similar picture 
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Figure 4.7: Smooth version of figure 4.6 with / = — 1 — tanhy 



2. what is happening in the smooth version is that the tangent becomes verti- 
cal, so |^ — > oo; 

3. if, as happens in some examples, z is the height of a fluid surface, then this 
is like a breaking wave but if z was e.g. pressure p(x, y) at a point in the 
plane, we would need a rule to help us choose which of the solutions is the 
physically correct one. We would then need to learn to live with solutions 
with discontinuities, interpreting them as "shocks" or "shock waves". 
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5 Classification of 2nd order linear PDEs 

In this section, we are interested in PDEs of the following form: 

a(x,y)u xx + 2b(x,y)u xy + c(x, y)u yy + f(x, y, u u x , Uy) . (5.1) 
principal part linear 

We know to call these second-order and linear, and we have seen the following 
examples in Mods: 

u xx + u yy = Laplace equation 
u X x — Uyy = wave equation if y — ct 
m„ — u y = heat equation if y = t/k. 

5.1 The idea: 

In this section, the key idea is to change coordinates so as to simplify the principal 
part. Thus 

(x,y) (ip(x,y),ip(x,y)); 

with 

d((p : ip) 

W^) = ^ " ^ ^ °' 

For the change in the partials, we calculate 

U X = U V if X + U^l/j X ; Uy = U^y + U^lfjy 

then 

u xx = u w Lp 2 x + 2u (p - l p(p x ip x + u^ipl + u^ xx + u^ xx 

U X y = U W tp x tPy + U^((p x tPy + 1p X (Py) + U^l]) x l])y + U^xy + U^ X y 
Uyy = U W (f 2 y + 2U lp ^(Py1py + U^lpy + U^ifyy + U^yy 

so that (5.1) becomes 

A u w + 2Bu lfyi p + Cu^p + F(ip, ip, u, u v Uip) = (5.2) 
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with 

A = atpl + 2b<p x (f y + cif 2 y "j 

B = a(fi X l/j X + b((p x lfjy + (Py^x) + CiPy^y > (5.3) 

C = mfe + 2b^ y + ci> 2 y J 
In a matrix notation (5.3) is (check!) 

(A B\ = fip x <f y \ fa b\ fip x ip x \ 

\B CJ \4) X ijyj \b Cj \(fiy ijyj 

so that 

(AC - B 2 ) = (ac - b 2 )(<p x 4> y - VWy) 2 - (5-4) 

(we could obtain (5.4) directly from (5.3) but the matrix notation makes it more 
perspicuous). Now (5.4) leads to a classification of second-order linear PDEs which 
is then invariant under change of independent variables: 

5.2 The Classification 

Second-order linear PDEs are classified into three types as follows: 

1. ac < b 2 hyperbolic: e.g. wave equation; 

2. ac > b 2 elliptic: e.g. Laplace equation; 

3. ac = b 2 parabolic: e.g. heat equation. 

We shall look at the classification in terms of the quadratic polynomial 

a(x, y)\ 2 - 2b(x, y)X + c(x, y) = 0. (5.5) 

Case 1: hyperbolic type 

So ac < b 2 and the quadratic has distinct real roots Ai, A 2 . 
We solve the first-order quasi-linear PDEs: 

(f x + Ai(x, y)ip y = 0; i\) x + A 2 (x, y)i\) y = (5.6) 

to find new coordinates if, ip; then from (5.3) and (5.5) A = C = 0; while from 
(5.4) B^O. Divide (5.2) by B to obtain the equation in the form 

+ G(ip, V>, UjU^u^) = 0. (5.7) 

This is the canonical form for a hyperbolic equation; ip, ip are characteristic vari- 
ables; curves on which ip or are constant are characteristic curves. We can often 
solve (5.7) explicitly. 
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5.3 Examples 

(a) 



^xx ^yy ^* 



We already know how to solve this, but let us apply the method. So 
a = 1, 6 = 0, c = -1, and A 2 - 1 = 0. 

We can take 

Ax = 1, A 2 = -1 

and solve (5.6) 
by 

(p = x — y ip = x + y. 
(There is clearly lots of choice at this stage.) The equation has become 

which we solve at once by 

« = /M + g(4>), 

a solution known from Mods, 
(b) An example with data: solve 

(x + y) 

XU XX - [X + y)U X y + yUyy + _ ^ (U a - My) = 



with 



u = -(x — l) 2 , u y = on y — 1- 



The quadratic (5.5) is 

xA 2 + (x + y)A + y = 
= (A + l)(xA + y); 

so choose 

Ai = -1 \ 2 = ~- 

x 

and solve 

ip x -ip y = 0; V* i> y = 
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by 

Calculate 
so that 



ip = x + y; ip = xy. 
u x = u lfi + yu^ 

v>xx — u<p<p ~\~ ^yUiptp ~\~ y Utptp 

^xy ^ifip ~t~ xu<pip ~\~ ytiipip ~t~ xytiipip ~\~ u^p 

^yy ^tptp ~t~ ^xu^ip ~\~ x Uipip 

(it isn't worth remembering the second-order chain rule, since it is readily 
obtained from the first-order one when needed.) 

Now the PDE becomes 

= x[u w + 2yu v4> + y 2 u^} 

-(x + y)[u w + (x + y)u^xyu^ + u^\ 
+y[u<p<p + 2xu vi , + x 2 u^} 
+ (x + y)u^ 
= (Axy -(x + y) 2 )u vi , 

so 

= 

and the solution is 

u = f(<p) + gty) = f(x + y)+ g(xy). 
To impose the data, calculate 

u y = f'(x + y) + xg'(xy) 

so on y = 1, 

u = f(x + l)+g(x) = -(x-1) 2 
u y = f'{x + l)+xg'{x)=0. 

Differentiate the first: 

f(x + 1) + g\x) = x - 1 
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and solve simultaneously with the second: 



g'(x) = -l, f'(x + l)=x. 

Integrate to find 

X ^ 

g(x) = -X + Ci f(x) = - -X + C 2 

and substitute back in u(x, 1): 

I( x _ l) 2 = - x + Cl + l( x + l) 2 _ ( x + i) + C2 

to find 

Ci + c 2 = 1. 

Finally 

m = + y "> 2 ~ ( x + y "> + 1 ~ xy - 

Case 2: elliptic type 

Now ac > b 2 so (5.5) has a complex conjugate pair of roots. Can we solve 

ip x + X(x, y)(f y = = i) x + A(x, y)^? 

Assume so, with solutions ip — (p, then A = C = 0, 5^0 and the equation 
becomes 

Introduce real and imaginary parts for tp as ip = ( + ir], (p = ( — irj to obtain the 
canonical form for an elliptic equation: 

u C(: + u vv + H((,r),u,u u v ) = 0, (5.8) 
which closely resembles the Laplace equation. 

Case 3: parabolic type 

Now ac = b 2 so (5.5) has a repeated root X(x, y). Solve tp x + X(x, y)<p y = for one 
new coordinate, and pick any ip with p x i J y—<P y i J x ^ as the other, then A = B = 
so, provided C ^ 0, we get the canonical form for a parabolic equation: 

+ G(u, tp, ^, «v) = °> 
which closely resembles the heat equation. 

59 



5.4 Example 

Classify and reduce to canonical form the equation 



x u xx + 2xyu xy + y u 



"yy 



0. 



The relevant quadratic is 

x 2 \ 2 - 2xyX + y 2 = = (xX - y) 2 
which has equal roots, so this equation is parabolic; A = | so solve 

<Px + -<P y = 0, by, e.g., <p=- 

X X 



and take, for example, ip — x. Calculate 

u. 



Uy 



1 

X 



so that 



u 



y q y i i <") y 



x;> 



x" 



u 



y 



xy 



X 6 X X 2 



u 



The equation becomes 
2y 



x 



r ~*y % 



+ 



+ 2xy 

so the canonical form is 



y 



X 3 W X 



1 

a; 2 



+ y z 



X 



2 tt'ipifi 



(5.9) 



x 2 u^ = (5.10) 



with general solution u = F(tp) +ipG(tp). In terms of the original variables this is: 

- F (;) +i °(;)- (5 - n) 

NB Very often, a question like this will be phrased in the form 'Classify and reduce 
to canonical form the equation (5.9) and show that the general solution can be 
written as (5.11)'. Therefore candidates for ip and ip are proposed by the question 
itself. 
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5.5 A warning example 



The type can change e.g. classify the equation 

^xx ~t~ U^yy 0. 

Then 

X 2 + y = 0, \ 2 = -y, 

and this is: 

• elliptic in y > 0, 

• parabolic at y — 0, 

• hyperbolic in y < 0. 

5.6 Type and data: well-posed problems 

We want to say something about the notion of well-posed-ness and its connection 
with type. Our examples are mostly based on knowledge acquired in Mods. 

A problem, consisting of a PDE with data, is said to be well-posed if the solution: 

• exists 

• is unique 

• depends continuously on the data. 

We can't be precise about what continuous dependence on the data means, but 
motivated by the general idea of continuity we take it to mean that a small change 
in the data (in some sense) leads to a small change in the solution (in possibly 
some other sense). So that if fi and f 2 are two sets of initial data and Ui and u 2 
the corresponding solutions then Ve 35 such 

1 1 A -/2II < 8 =>• IK - W2II < e- 

To be precise we need to specify what 'close' means, i.e. to specify a notion of 
distance between two functions |K — u 2 \\ and this takes us beyond this course. 
Examples might be the supremum of K — u 2 \ over some domain or f K — u 2 \ 2 
but there are many more and that for u will in general be different to that for /. 

We look at some examples from Mods: 
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(a) The IVP and IBVP (initial-boundary-value problem) for the wave equation 

u xx -u yy = (ct = y). 
For the IVP, we know the solution is 

u = f(x + y) + g(x-y) 

= 2 + y) + v( x - v)\ + 2 j i>( s ) ds , 

where the data are u(x, 0) = ip(x) and u y (x, 0) = ip(x). This is d'Alembert's 
solution of the IVP: it exists, and is unique and, intuitively at least, a small 
change in ip, tp gives a small change in u (in fact one can be a bit more precise 
here). So this problem is well-posed. 

For the IBVP consider the set-up: 

u(x, 0) = f(x), u y (x, 0) = g(x) < x < L 

u(0,y) = = u(L,y). 
So the boundaries are at x — 0, L. This IBVP is solved with Fourier series 

Ennx ( rnry nny\ 

sin — — [a n cos — — + b n sin — — J 



n 



with a n , b n uniquely fixed by /, g. Again this problem is well-posed (with 
an appeal to intuition for continuous dependence on the data). 

(b) The BVP for the Laplace equation 

Do this first with data at the sides of a square, so < x, y < a with 
w(0, y) = u(a, y) = u(x, 0) = 0; u(x, a) = f(x). 

Consider separable solutions u n = sin ^ sinh ^p, then 

Enirx niiy 
a n sin smh 
n n 

and 



u 

n 



f(x) = a n sinh(n7r) sin 
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nnx 



which determines the solution as a Fourier series. 

Now a different BVP, with data at the circumference of the unit circle: 

on r = 1, u — f(9) 

and in polars 

1 1 

U„. + -U r + —zUqq = 0. 



The separable solutions are 



(Ar n + (C cos n6 + D sin n0) 
A + B log r, n = 

Regularity at r = implies 

j X 

1 

and the boundary value at r = 1 requires 

-ao + ^^( a n cosnO + b n sin n6>) = f(9) 
which again is solved by Fourier methods. 

For these two cases at least, the BVP for the Laplace equation is well-posed 
(uniqueness was shown in Mods). It is plausible, but beyond our scope, to 
show this in general. 

(c) The IBVP for the heat equation 



on the semi-infinite strip where y = t > and < x < L, and data 
u(x,0) = f(x); u(0,y) = = u(L,y). 



The relevant separable solutions are 



nirx n 2 vr 2 t 



sm — — e l 



— 



L 

so that 



u = a„ sin - - -< 



U7TX n it 

L 



2_2 

t 



n 
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and the initial value requires 



/(*) = E 



nirx 



a n sin 



L 



which is solved by Fourier methods. The solution exists provided the series 
for u converges which it will do for positive t. However, note that for negative 
t the exponentials grow rapidly with n and there is no reason to expect 
existence. Uniqueness for this problem was done in Mods, and we can argue 
plausibly for continuous dependence on data, but we only get well-posed-ness 
forward in time. 

(d) What then is not well-posed? We give a few examples: 

• BVPs for hyperbolic 

e.g. u xx — u yy = on the unit square with data 



Recall this data gave a well-posed problem for the Laplace equation, 
but here there is no solution at all if / ^ (try the Fourier series) while 
if / = 0, then sin no; sin ra/ will do, for any n. 

• IBVPs for elliptic 

e.g. u xx + u yy = on the semi-infinite strip < x < 1, y > 0, with data 



This data gives a well-posed problem for the wave equation. If we try 
for separable solutions here, we have u n = sin nirx cosh nny so 



u(0, y) = u(l,y) = u(x, 0) = 0; u(x, 1) = f(x). 



u (0, y) = u(l, y) = 0, u(x, 0) = 1, u y (x, 0) = 0. 



u = 




n 



Initial conditions need 




whence 



a n = n even 
4 

= — n odd, 
nn 



and then 




n 
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which does not converge for any y > (because the cosh terms grow 
rapidly with n) - there is no solution (strictly speaking, we've only 
shown that there is no solution of the form considered; we need more). 

• the BVP for the heat equation is not well-posed, but we won't show 
that. It's easier to see that going the wrong direction in time is not 
well-posed: consider the particular solution 

1 

u(x,y) = — e 4 y . (5-12) 

Vy 

Start with this solution at any positive value of y, say at y — e, then 
it is analytic in x and its graph looks like the normal distribution; now 
evolve it towards negative y; at any nonzero value of x it tends to zero 
as y tends to zero, but at x — it diverges, and it makes no sense at 
all for negative y. The solution has failed to exist within a 'time' e. 

Again, it is beyond our scope to prove it in this course, but these different 
behaviours are universal for the different types of second-order, linear PDEs. 
In tabulated form, which problems are well-posed? 





IVP 


IB VP 


BVP 


Hyperbolic 


yes 


yes 


no 


Elliptic 


no 


no 


yes 


Parabolic 


yes 


yes 


no 



where the 'yesses' for parabolic are forward in time only. 



5.7 Bessel Functions and Legendre Polynomials 

We end this chapter with two odd topics on separation of variables, which don't 
fit anywhere else. 

5.7.1 The vibrations of a circular drum in polar coordinates 

Suppose the undisturbed drum occupies the disc r < a, z = in plane polar 
coordinates, and the vertical displacement of the disturbed drum head is given by 
z = u(r, 9, t). Then u satisfies the wave equation in the form 

1 n2 11 

— u tt =Vu:=u rr + -u r + -^uee- 
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We seek separable solutions 



u = R(r) 6(0) T(t) 

with 

T = C ° S cot 
sm 

so 

c z r r z 

where prime on a function of one variable means derivative with respect to that 
variable. Divide by RQ and separate to find 

, ,'R" 1 Rl co 2 \ 6" , £ 

r ( — H — H — 7 = — — therefore = const. 

R r R & ) 6 

but we want 6 periodic so the separation constant must be n 2 for integer n: 

, C ° S n6 n > 1 
9 = < sm — 

1 n = 

1 / /w 2 n 2 



r \ c 2 r 2 



r ■ -/; + i — - — ) /; o. 

Redefine p — -r to obtain 



d2 « + If + f 1 _4' |fl = . (5.13) 



<ip 2 p dp \ p 2 

This (second-order, linear ODE) is Bessel's equation of order n. We may seek a 
series solution of the form 

oo 

R = p «Y^a k p k 

o 

in terms of constants a and a k , by substituting this into (5.13). This is straight- 
forward to do, but not part of this course. The outcome is that, up to scale (since 
any constant multiple of the solution has the same property), there is a unique 
solution which is finite at p = 0. With a conventional choice of the scale, this 
solution is called J n (p), the Bessel function of order n. The series begins: 



p n 



1 " 71 ^ T\ + °(P 4 
4(n + l) Kl 
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(it's easy to check by substituting into (5.13) that there is a series solution like 
this). 

Given the series, one can plot the function to obtain graphs like figure 5.1. 




Figure 5.1: The first two Bessel functions 



For the problem of the drum, we want u — at the rim r = a, so R(a) = and 
so J n (^) = 0. Thus ^ = a nm , where this is the m-th zero of J n , which can 
be found numerically from the series solution. Therefore the allowed frequencies 
of oscillation of the circular drum are u nm = -a nm . The disturbance u is a 
combination of terms like 




5.7.2 Laplace's equation in spherical polars 

Recall spherical polars are related to Cartesian coordinates by 

x = r sin 9 cos tp 
y = r sin 9 sin tp 
z = r cos 9. 

and then Laplace's equation becomes 

2 1 1 

V 2 m := u„ + -u r + - , n (sm9u e )e + . 2 u w = 0. 
r r 2 sm6 r 2 sm 9 

We could separate this in general as u(r, 9, tp) = R(r)Q(0)Q(<p), but for simplicity 
we shall suppose u v = 0, so that u = R(r)Q(9). Substitute into Laplace's equation: 

q( R " + * a) + ^ J-(sin0e')' = o, 

\ r J r 2 sm9 
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where prime on a function of one variable again means derivative with respect to 
that variable. Divide by RQ and separate to find 

r ^-(R" + -R!) = -t^— - (sm9ey = const. = -A 
R K r J Gsintf v ; 



so that 



R" + -R' + \r = (5.14) 



(sin 00')' - A sin #6 = 0. (5.15) 
To simplify (5.15), put /i = cos#, so that sin#J^ = —(1 — /i 2 )^ and 

iH (i -" 2) fH e - (516) 

which is Legendre's equation. As with Bessel's equation, we seek series solutions 
= ^afc/i fc with constants a^. This time we find that the series converges for 
\x = ±1 (which it must for solutions to be defined at 9 — 0, n) if and only if it 
terminates, i.e. is a polynomial. If the polynomial has degree n, then necessarily 
A = — n(n + 1). These solutions are Legendre polynomials P n and the first few can 
be taken to be 

Po = l| Pi = K P2 = \(Sfi 2 -l), .... (5.17) 



For (5.14), this leaves 

R"+'- 

for which solutions are r n , and r~ n_1 , so the general solution for u is 



R" + 2 -R> - n { ^AR = 



u 



= E(^ n + ^r) p «( C0 ^) ( 5 - 18 ) 



An example: 

Consider a BVP for the Laplace equation on the region between concentric spheres 
of radius a and b, with < a < b. In spherical polars suppose the solution is u(r, 9), 
for the data 

u(a, 9) = 1; u(b, 9) — b cos 9. 

The solution is unique by Mods work, and we expect well-posed-ness. From (5.17) 
and (5.18) we expect a series with just the first two Legendre polynomials Pq and 
Pi- 

u = A + — + \A x r + -4 J cos#, 
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and then the boundary conditions require 



A) + — + ( A ± a + ^ ) cos9 
a V a 1 



For these to hold at all 9, we obtain four equations in four unknowns: 

Ao + lt = 1 

+ = 

a 2 

Aib + ^ = b 
b 2 

from which a unique solution follows by elementary algebra. 
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6 Laplace and Fourier Transforms 



In this section, we introduce a new idea, the integral transform which can convert 
differential equations into algebraic ones (which are typically easier to solve). We 
shall consider two such transforms. 

6.1 The Laplace Transform 

of a function f(t) is 

PCX) 

f(p) = / f{t)e~ pt dt. 
Jo 

Note that 

• the over-bar does not indicate complex conjugate; 

• a convenient notation is to write / = £[/]; 

• we shall often allow p to be complex and then be interested in where / is 
defined in C; 

• if 1/(01 ^ Me ct for some positive M, c then the integral is defined at least 
for some p (in fact for all complex p with Re(p) > c); such / is said to have 
exponential growth and examples of such / are any polynomial in t, or any 
exponential but not e.g. e* . 

6.2 Examples: a standard list 

Check the following standard list of examples of Laplace transforms, with corre- 
sponding domain of definition: 





/(*) 


Kp) 


Domain of definition 


1. 


1 


i 


Re(p) > 


2. 


t 


I 

p2 


Re(p) > 


3. 


t n 


n! 

pn+l 


Re(p) > 


4. 


e at , aeC 


p—a 


Re(p) > Re(a) 


5. 


sinwt 


U) 

P 2 +U! 2 


Re(p) > 


6. 


cos out 


V 

p 2 +ui 2 


Re(p) > 


7. 


sinh at 


a 

p 2 —a 2 


Re(p) > \a\ 


8. 


cosh at 


P 

p 2 —a 2 


Re(p) > \a\ 
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(to save work, note that (5) - (8) follow from (4).) 



6.3 Properties of the Laplace transform 

1. Linearity: C[\f + jig] = A£[/] + fJ>£\g], for A,// G C; 

2. If g(t) = e at f(t), then g(p) = f(p - a); 

3. If g(t) = /'(f), then g(p) = pf(p) - /(0); 

4. Ug(t) = fif(8)d8, then g(p) = M ; 

5. Define the convolution f * # of / and g by / * <?(f) = J Q * /(f — s)g(s)ds, then 
£[/ * </] = £[/]£[<?] or = / g- 

6. If =«/(*)» then ^(p) = -^/(p). 

Proofs 

1. easy; 

2. $(p) = / °° g(t)e-*dt = / °° f(t)e at e-*dt; 
= /o°°/(Oe- (p - o) ** = /(p-a); 

3- 3(p) = JZ°f(t)e-*dt = \f(t)e- pt \?+Pj?f(t)e-*dt 

= -/(0)+p/(p); 

4. 4= (3.) 

5. £[/ * «?] = / °° £ /(f - S )^(*)dse-**dt 

now interchange the order of integration: 
= J? F f(t - 8)g( 8 )e-*dtd8 

change variable, eliminating t in favour of w = t — s 

= Jo°° Jo°° f(u)g(s)e-P s e-P u duds 

= f(p)g(p) = C\f\C\g]. 
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6- if(p) = if O °f{t)e->*dt 

= -f CO tf(t)e-*dt = -g(p)- 

Second derivative: corollary to property 3 

If g{t) = fit) then g(p) = p 2 /(p) -p/(0) - /' (0). 

(Exercise: introduce h — /', so g — h! and use item 3.) 
Now let's see the Laplace transform in action: 

6.4 Examples 

(a) Solve the IVP 

x" - 3x' + 2x = 4e 2 *; x(0) = -3; x'(0) = 5 

for x(t). 

Do C to the ODE and use property 3 and its corollary: 

4 



(p 2 x + 3p - 5) - 3(px + 3) + 2x 



p-2 

where the RHS comes from our standard list; so 

/ 2 « 4 n -3p 2 + 20p - 24 

x p 2 - 3p + 2 = - 3p + 14 = — ^ £ 

p — 2 p — 2 

which can be solved for x as 

-3p 2 + 20p - 24 
X ~ (p-l)(p-2)2 • 

Split into partial fractions 

-7 4 4 
+ t + 



p-1 p-2 (p-2) 2 ' 

and use the standard list to give 

x = -7e* + 4e 2t + 4te 2t 

where the last term needs C[e 2t ] = so = ~~ ^^[ e2< ] = £[ te2< ] by 

property 6. 
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(b) Solve 

x " + uj 2 x = f(t) x(0) = a, x'(0) = b 
(this is the equation for a harmonic oscillator with a driving term f). 

Laplace transform of both sides gives 

(p 2 x — ap — b) + cu 2 x = f 

so 

x{p 2 + w 2 ) = f + ap + b 

and 

- - "'' I b , / 



G» 2 + w 2 ) (p 2 + w 2 ) (p 2 + w 2 )' 
The last term, by property 5, is a convolution, so 

x = a cos H sin + / * sin cut 

u 

or 

b If 1 
= a cos ut H — sinwtH — / /(s) sinw(t — s)ds 
u to J 

(we could have solved this by the Green's function method, and it's worth 
comparing the answers.) 



(c) an ODE with non-constant coefficients (this doesn't always work). 

tx" + 2x' + tx = x(0) = 1, 
so from the ODE, for a regular solution we also need x'(0) = 0. 

Laplace transform the equation: 

c[tx"\ = -^£[A = -^(p 2 x-p) 

S ° d d 

— —(p 2 x — p) + 2(px — 1) — —x = 
dp dp 

and 

- ( i + *>)|-i = o 
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so that 



dx 1 



£[sin t] 



dp p 2 + 1 



from the standard list. We've reduced a second-order ODE in t to a first- 
order one in p but we can solve this using property 6, to find tx = sint or 



In these examples, we have needed some luck to get back from x to x, i.e. to invert 
the Laplace transform. How can this be done in general? We shall postpone this 
question until after the discussion of the second integral transform. 

6.5 The Fourier Transform 

Much of the study of this transform can proceed by analogy with the previous. 
First the definition of the Fourier transform of f(t): 



x = 7 sint. 




Then note 



• existence needs only \f(t)\dt < oo (this is a fact from the Integration 



course) 




• some authors include a factor -4= in the definition; 



• we shall write / = F[f]. 



6.6 Properties of the Fourier Transform 



1. linearity: 



F[\f + fig] = XF\f\ + nF[g] 



2. if g{t) = fit), then g{u) = iuf(oj) (so no "/(0)" as compared to the Laplace 
transform) ; 
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3. 




lution agrees with the earlier one if / = g — for t < 0, but otherwise is 
different.) 

Proofs 

Exercises. 

6.7 Inversion of the Fourier Tranform 

This is given by the Inversion Formula: 



where the expression on the left is the average of the limits of / approaching t 
from above and below. If / is continuous at t this is just f(t) . 



This will be just a sketch. I shall indicate where input from Complex Analysis or 
Integration is required. Recall from Complex Analysis that 




(6.1) 



Proof 




then 





— oo 



-R 




(OK to change the order of integration?) 




(OK at t = s?) 




75 



= 2(/ 1 + / 2 ) 

where 

/"* sin R(t — s) , . . , /"°° sin if . . u . , 

'> = L—j^r-H'V" I — ^-«>*' 

with t — s = -|, and 

r°°smR(t-s) £l , , psinw^ 

J 2 = / >-f(s)ds= / f(t + -)du 

Jt t-s J u R 

with s — t = -|, so 

£«■"/<«**» = 2 f ^ (/('- 1) + /<« + 1>) 

Now take i? — >■ oo and use the integral with which we began (is the limit OK?). 
QED. 



6.8 An example 

Invert 



/= 2 



1 + w 2 

The Inversion Formula (6.1) gives 



fit) = — lim / -duj 

J w 2vr R^oo J_ R 1 + u 2 

which we shall evaluate by closing the contour with a semi-circle of radius R 
centered at the origin. If t > we choose the semi-circle in the upper-half-plane, 
and conversely if t < we take it in the lower-half-plane. Call these contours 
C + and C- respectively. Poles of the integrand are only at ±i so for positive t, 
Cauchy's integral formula gives: 

e iujt 1 

rdoj = 2iri x — e~* = 7re~*, 



c+ 1 + to 2 2i 
while if t < 0, we obtain 

1 



/ = -2ni x {-i-e 1 ) = rre 1 
JC- 2« 



the extra minus sign coming from the fact that the contour is traversed clockwise. 
Putting these together we obtain f(t) = e~'*L 



76 



6.9 Inversion of Laplace Transform 



From (6.1), we can obtain an inversion formula for the Laplace transform as follows: 

1 1 p+iR 

-(f(t-) + f(t + )) = —hm / e*f(p)dp (6.2) 

where 7 is chosen as follows: if \f(t)\ < Me ct then 7 > c. This means that the line 
along which the integration is carried out is to the right of any singularities of /. 
As before, if / is continuous at t then (6.2) gives f(t). 



Proof 



This uses Fourier transform and (6.1). Given /, define 

9(t) 

Then 



e~^f(t) t > 
t < 



poo 

g(u)= / e--*e-^/(*) = /(7 + H 
Jo 

and, by (6.1) 



2 

so 



i(0(t_)+0(t + )) = J- lim / e^f( 1 + iuj)dw 

I In R-^co J_ R 

\{f{t_) + /(*+)) = ±- lim f R e (TH»>*/( 7 + iu)dw 

I Z7T R^oo J _ R 

■y rj+iR 

= lim / e pt f(p)dp where p = 7 + iu, 

2-Kl R->oo J iR 



>7-iR 

as required. QED 

When it comes to examples, we often evaluate the integral by closing the contour 
to the left. 



6.10 Examples 

(a) Invert 



ft?) = 2i 1 TV 
p 2 (p- 1) 

77 



This example can be done by partial fractions and the standard list, but we 
shall use the Inversion Formula. Poles of / are at and 1 so we need 7 > 1, 
and we shall close the contour to the left, with an arc of a circle centered at 
the origin. Write I\ for the straight part of the contour, T 2 for the curved 
arc and T for the union. By (6.2), at points of continuity of /, 

fit) = l]i m / e pt J V - 
and we claim that, as R — > 00 

2ni x sum of residues. 

Res\ = -(l + t), 
-(l + *) + e*. 

(b) Invert f(p) = ^. 

This has a branch point at the origin. We proceed as before, closing the 
contour to the left with a circular arc centered at the origin, but we need a 
'key-hole', excluding the negative real axis and the origin. 




so consider 

oPt d _l_ 



r p 2 (p-l) 
Calculate 

Res\i = e*; 

and take the limit to find that / = 
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Then we may define p 1//2 = r l / 2 e l6 ^ 2 for —n < 9 < it and consider the integral 
fr ePt ~T72- With points as labelled on the diagram and DEF a circle of radius 
e, we claim that 




To evaluate this integral, put rt = s 2 , 

1 r° _,2sdst 1 / 2 
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= — -= / e ds 

TTyt J-oo 

which is a standard integral: 

1 

"~ Vtn' 

With the aid of this example and property 6 of the Laplace transform, we 
can also invert p~ 3//2 , p~ 5//2 , etc. 

The next set of examples show how transforms can help with PDEs as well 
as ODEs. 

(c) Consider the following IBVP for the heat equation 
in x > 0,t > 0: 

u(x, 0) = initial condition 

—ku x (0,t) = Q 
u — > as x — ?■ oo 
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I boundary conditions 



Find u(0,f) for t > 0. 



The problem corresponds to a semi-infinite bar, initially at zero temperature, 
which has a constant rate of heat-flow into it at the origin. Heat will flow in 
and pass along the bar and in particular the end at the origin will heat up. 
We are asked for the temperature of this end. 

We Laplace transform in t only (we only care about t > 0): 



u 



poo 

(x,p) = / u(x,t)e~ pt dt. 
Jo 



Then 



£[u t ] = pu — u(x, 0) = pu 



roo 

Jo 

This is an ODE in x for u which we can solve: 

u = A(p)e x V^ + B(p)e- x V^ 
but u — > as x — )• oo. Therefore, A = and 



u x {0,p) 



Q 
k 

Q 
k 



u x (0,t)e" pt dt 



- pt dt 



--e 
p 

Q_ 

kp 



- P t 



-\I- K B(P) 



so 



whence 



B = 



QV^ -3/2 
k P ' 



u(x,p) 



QV^ p -3/2 c -xy/^Jk 

k 
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Inverting this would give us u(x, t), but we are only asked for u(0, t). For 
this 



k 



To invert this, go back to example (b): 

1 



C 



1 

y/P' 



so 



d 
dp 



(P- 1/2 ) = -^ 3/2 



1 

2 1 
ds 



= -C 



using property 4 of Laplace transform, 



-C 



performing the integral, so that p 3//2 inverts to 4< / -, and 



u{0,t) 



AQ rid 



k \ 7Y 



(d) An IVP for the heat equation in an infinite bar: 

u xx = u t ] —oo < x < oo, t > 

u(x,0) = f(x). 
This time we Fourier transform in x: 



/oo 
u(x, t) 
-oo 



e~ lwx dx 



so 



an ODE we can solve: 



but 



u{u,t) = A(w)e" w * 
■u(cu, 0) = /(w) therefore A — f 
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and 

u(u,t) = /He"" 2 ' = f(u)k(w,t). 
This is a product of Fourier transforms, so corresponds to a convolution: 



/oo 
K(x-y,t) f(y)dy (6.3) 
-oo 



where 

AT(u;,t) = e- wH . 

By the inversion formula (6.2) 



1 

iT(x,t) = — lim / 

2,71 R—>oo J_ R 

so with s = tjj\ft this is 



| pH 

— lim / e- s2+isx/Vl ds. 



We can evaluate this integral by the following Lemma. 
Lemma: for real a 



/oo poo 
-oo J —oo 



Proof 



By Cauchy's integral formula 

je- {z - ia)2 dz = 

where T is the rectangle with vertices at (R, R+ia, —R+ia, —R) for positive, 
real R. Now take R to infinity and show that the integrals along the short 
sides tend to zero, so that 



/R pR 
e- {x ~ ta)2 dx = lim / e~ x2 dx = ^ 
-R R-^too J 



-R r R 

R-^c 

With the aid of this Lemma, we have 



i -x 2 /At 

K(x, t) = ■= ■ ^e~ x2 ' u = (6.4) 
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and then 



1 r°° 

u(x,t) = —= e-^-y) 2 f{y)dy. (6.5) 



It's worth making a few points about this calculation: 

• note the resemblance of (6.3) to Green's function expressions from sec- 
tion 2, like (2.12) or (2.16); 

• note that K(x,t) in (6.4) has appeared before (apart from a constant 
factor), in (5.12); it is a solution of the heat equation which, at t — 0, 
is zero everywhere except at the origin, where it is infinite; from its 
resemblance to the normal distribution, we know that the integral over 
x of K is one at all times - intuitively it is like a unit point source of 
heat at time zero spreading out; 

• we can think of (6.5) as expressing u(x,t) as a superposition of point 
sources distributed according to the initial temperature given by / and 
spreading out; 

• as we saw in the discussion of well-posed-ness, the solution is defined 
for positive t only. 

(e) The Poisson kernel 

Something similar to the last example is possible for the BVP for the Laplace 
equation. Consider the problem: 

u xx + u yy = , — oo < x < oo , y > 

u(x,0) = f(x) 
u — > as y — > oo. 

So this is the BVP with data on the real axis. We Fourier transform with 
respect to x: 



/oo 
e~ luJx u(x, y)dx 
-oo 



and then use the Laplace equation and property 2 of the Fourier transfor- 
mation to find 

Uyy = u 2 u, 

with u(uu, 0) = f(oj), and u — > as y — > oo. The solution of this ODE for u 
is 

u = A(cu)e~^ y + BeM", 
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taking account of the fact that u may be positive or negative, and then the 
boundary conditions in y imply that 



u(u,y) = f(u)e~^y. 



This is a product, so the answer will be a convolution. For the inversion of 
the second factor, we consider the integral 



1 



2tt 



R 



— / e^-^diu 



R 



2tt 



— / e^+^du + — / e^ ix ~ v) dw 



R 



1 

2^ 



e u)(ix+y) 

ix + y 



-R 



2tt 



1 



R 



e ui(ix-y) 

ix — y 



R 



and for large R: 



1 



2n \ix + y ix — y J n(x 2 + y 2 ) 



■= K(x,y). 



Therefore 



/oo 
K(x - s, y)f(s)ds 
-oo 

m 



(x — s) 2 + y 2 



ds. 



(6.6) 

(6.7) 
(6.8) 



The function K is called the Poisson kernel. 

The remarks after the previous example carry over here: (6.7) is like a Green's 
function formula; K in (6.6) is a solution of the Laplace equation singular at 
one point; (6.8) expresses u as a superposition of these elementary solutions 
along the boundary. 
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